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Abstract. We apply transition state theory to coupled Gaussian wave packets and 

calculate thermal decay rates of Bosc-Einstcin condensates with additional long-range 
interaction. The ground state of such a condensate is metastable if the contact 
interaction is attractive and a sufficient thermal excitation may lead to its collapse. 
The use of transition state theory is made possible by describing the condensate within 
a variational framework and locally mapping the variational parameters to classical 
phase space as has been demonstrated in the preceding paper [A. Junginger, J. Main, 
and G. Wunner, preceding paper, submitted to J. Phys. A]. We apply this procedure 
to Gaussian wave packets and present results for condensates with monopolar 1/r- 
interaction comparing decay rates obtained by using different numbers of coupled 
Gaussian trial wave functions as well as different normal form orders. 
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1. Introduction 

Since their first experimental realization in 1995 yj Bose-Einstein condensates (BECs) 
have become an active field of theoretical and experimental investigations. Moreover, 
BECs with additional long-range interaction are of special interest, because the 
interactions can be tuned from predominantly short-range to the dominance of the long- 
range interaction by manipulating the s-wave scattering length via Feshbach resonances. 
The latter allow varying the contact interaction in strength as well as in sign. In case of 
a negative scattering length, i.e. an attractive interaction, the ground state of the BEC 
is metastable so that the condensate may decay by collapsing after a sufficient thermal 
excitation. 

The thermal decay rates of BECs without long-range interaction have already been 
estimated by Huepe et al [2|[3] within a simple variational ansatz of a single Gaussian 
wave function. This approach can, of course, also be applied to condensates with long- 
range interaction, but it will, because of its simplicity, only yield qualitative results. 

In this paper, we present an improvement by using an extended variational ansatz 
with coupled Gaussian trial wave functions, which are described in the framework 
of a time-dependent variational principle and which have proven their capability to 
reproduce the numerically exact results or even to exceed them (4|[5]. Within this 
variational ansatz, the BEC exhibits two stationary states, one of which corresponds 
to its metastable ground-state, while the other one is an excited state of saddle-centre- 
. . . -centre type. The decay rate of the condensate can, thus, be calculated by means of 
transition state theory (TST) [6] because the collapsing BEC has to cross this saddle in 
the subspace of the variational parameters and the decay rate is given by the fiux over 
the saddle. 

Although classical TST requires the knowledge of a Hamilton function H{q,p) given 
in phase space coordinates q,p, its application is made possible by locally mapping the 
variational parameters to action variables of the classical phase space in the vicinity 
of the fixed points. This local mapping to phase space is performed with the use of a 
normal form expansion of the equations of motion determining the time evolution of the 
variational parameters as well as the respective mean-field energy functional. 

For systems with known classical Hamiltonian this procedure has been shown to 
reproduce the decay rates of the classical and the quantum normal forms in the limits 
of narrow and broad wave functions, respectively |7], and, moreover, it well applies to 
systems where such a Hamilton function is not directly accessible as it is the case for 
the variational approach to BECs with coupled Gaussian wave functions. 

In order to demonstrate the applicability to BECs and to calculate their decay 
rates, our paper is organized as follows: First, we review the description of the 
variational ansatz in the framework of a time-dependent variational principle as well 
as the procedure of locally mapping the variational parameters to phase space. Then, 
we illustrate the calculation of the thermal decay rates by applying TST and at the end 
present and discuss the results for BECs with monopolar 1 / r-interaction. 



TST for wave packet dynamics. II. Thermal decay of BECs with long-range interactions 



2. Theory 

We consider a condensate consisting of bosons which exhibit an additional long-range 
1/r-interaction as has been proposed by O'Dell et al [s]. Such systems have not yet 
been experimentally realized but because of the spherical symmetry of this interaction, 
they form an important model system. 

At ultra-low temperatures this quantum gas can be described by a single wave 
function ip{r, t) whose time evolution is given by the extended Gross-Pitaevskii equation 
(GPE) 

Hi;ir,t) = idMr,t), (1) 
where the mean-field Hamiltonian 

H = -A + Vt + V, + V^ (2) 

describes the interaction with an external trapping potential Vt, the inter-atomic contact 
interaction Vc with the s-wave scattering length age, and the long-range monopolar 
interaction Vm- 

V, =NS'r^ (3) 
K =87^a,,|V'(r)|^ (4) 

V^.--2/dV^. (5) 



To obtain the GPE in this dimensionless form, we introduce "atomic" units [9] with 
the help of the constant u which determines the strength of the attractive monopolar 
inter-atomic interaction V{r,r') = —u/\r — r'\ Lengths are measured by means 
of the "Bohr radius" = h?/{mu), energies in units of the "Rydberg energy" 
Eu = ^^/(2"^a^), and times in units of tu = h/Eu- In addition, we apply a particle 
number scaling according to 

r ^N-^a^r, (6) 

(7) 

t N-hut, (8) 
E ^N^E^E, (9) 

^ (a„iV)-=^/V, (10) 
which eliminates the explicit occurrence of the particle number in the GPE. 

2.1. Application of transition state theory to the Gross-Pitaevskii equation 

The usual way of solving the GPE, equation ([T|, is either by performing an imaginary 
time evolution on a grid or by integrating it outward with initial values for the wave 



function and its derivative 10 . The former method can also be applied to more 



general geometries and more compUcated interaction potentials such as the dipole-dipole 
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interaction, while the latter is limited to the case of effectively one- dimensional systems 
like the radially symmetrical BEC investigated in this paper. 

In order to apply TST it is furthermore essential to precisely define the transition 
state of the system and therefore to find the unstable stationary solution of the GPE. 
However, since an imaginary time evolution will only yield the ground state of the system 
this method cannot be applied. Thus, the approach may be limited to the simplest case 
of a monopolar BEC in a radially symmetrical trap. 

If the stable (s) and the unstable {u) solution of the GPE are found the 
corresponding local normal form Hamiltonians are required in addition. Their lowest- 
order quadratic approximations will take the form [6] 



H^^\q,p)=Hi'^ + Y,-^uf^V,q,. (11) 

d 

H(-\q,p) = + Xp,q, + J^ic^f (12) 

respectively, where A is the eigenvalue corresponding to the decay channel of the BEC. 
Assuming their knowledge, the evaluation of the respective phase space integrals 11 is 
trivial and one obtains the decay rate 

r = '^e'"^''^ ) n (13) 

i=2 

which is a generalization of the formula given in reference [s] to c? degrees of freedom. 

The frequencies cjj*'"'' can, in principle, be obtained by the Bogoliubov-de Gennes 
(BdG) equations. However, this would present difficulties since they exhibit an 



unbounded spectrum so that the decay rate given by equation (13) will, in general, 
either diverge or vanish. Moreover it is questionable, to what extent a local quadratic 
approximation will yield appropriate results after integration over the whole phase space. 



An extension of the Hamiltonians (11) and (12) to higher-order terms would 



require the treatment of small perturbations to the solution of the GPE in higher-order 
approximations than the linear one resulting in the BdG equations. In conclusion, a 
numerically exact approach for the application of TST to the GPE does not appear 
promising. 

All these problems can, however, be circumvented if one treats the GPE within a 
variational framework: This is, in principle, not limited by any restrictions concerning 
the geometry of the system and the number of degrees of freedom, respectively. 
Moreover, it can be applied to various interaction potentials, allows a rather simple 
determination of the stable and unstable stationary solutions of the GPE and with it a 
precise definition of the transition state. Even the problem of unbounded spectra does 
not occur because one makes use of a finite set of variational parameters. Furthermore, 
an extension to higher-order normal form Hamiltonians can be carried out since the 
respective terms are related to numerically easily computable higher derivatives of the 
equations of motion which determine the time evolution of the variational parameters. 
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2.2. Variational approach to the Gross- Pitaevskii equation 

In the following, we will describe the condensate by means of a variational ansatz 

2P{r, t) = ij{r, z{t)) = exp(Ar' + 7i) (14) 

1=1 

in the form of a radially symmetrical Gaussian wave packet. Here, the complex and 
time-dependent variational parameters Ai and ji determine the width, the phase and 
the weight of each Gaussian, and we summarize all these parameters in the vector 
z = (zi, . . . , 27vJ^ with Zi = (Ai, 7i). 

An approximate solution of the time-dependent GPE ([T]) within the parameter 
subspace of the wave function ( 14 ) is given by the McLachlan variational principle [l2] 

I = \\i(f) - HtpW = mm. (15) 

where the quantity I is minimized with respect to and (p = ip is set afterwards. Its 
application to the parametrized wave function ( 14 ) yields the set of first-order differential 
equations 



13 



Kz = -ih, (16) 

which determines the time evolution of the variational parameters, and where the matrix 
K and the vector h are defined by 

hm = [d'r (P^Yh^. (18) 



The mean-field energy of the condensate is given by the expectation value 

E{z) = j dV r{r) (^-A + K + ^(K + Kn)) ^{r) (19) 

and is also a function of the variational parameters. 

For the ansatz with coupled Gaussian wave functions (14), all integrals occurring 
in equations (17)-(19) can be calculated analytically. However, since these calculations 
have been illustrated in detail elsewhere and are not subject of this paper, we refer the 
reader to references lillsl for their evaluation. 



2.3. Mapping to phase space 

The application of TST in phase space [6] requires knowledge of a (local) Hamilton 
function which describes the dynamics of the BEG in the vicinity of the unstable fixed 
point corresponding to the "activated complex". Such a Hamiltonian can easily be 
obtained even globally if one uses a single Gaussian to approximate the BEG's wave 
function [9]. In contrast to that, this is not possible in the case of coupled wave functions, 
as we use them in this paper. We therefore apply the procedure presented in reference IT] 
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to construct a local Hamilton function in the vicinity of the fixed points and in this 
section give a short overview on the steps performed (see reference [T] for details). 

To obtain the local Hamilton function, which, equivalently to equations (16) and 



(19), describes the dynamics and the energy of the system, one first Taylor expands the 
equations of motion (16) in the vicinity of a fixed point Zq up to the order rzmax and 
splits the expansion into its real and imaginary part. This yields the real vector field 

^ arXx) (20) 



X 



a[x] 



n=l 



with X = {Re{zi — zqi), lm{zi — zqi), . . . , Re^za — zod), ^^{zd — zod)) being the deviation 
of the variational parameters from the fixed point and a„(a3) summarizing all terms 
homogeneous of degree n. 

In the next step, equation (20) is diagonalized with respect to its linear part 
ai{x) = Ai ■ X and to further "simplify" the higher-order terms, a near-identity 
transformation x y given by (cf. reference (l4]) 

X = (l),{y), X = 0,=o(2/) = y (21) 

is performed. Here x = (f)^{y), which gives the identity transformation for e = 0, is a 
solution of the differential equation dx/de = g{x) and for an appropriate choice of the 
generating function g{x) brings equation ( [20| ) into the form {i = 1, . . . , d) 

E/5-(2.-i)^SiT^^r-' ll('^2,-iX2,r^, (22) 



X2i 



•m(2i) 



X. 



2i-l 



1X2 j 



(23) 



This form is due to the fact that the eigenvalues of equation (16) occur pairwise with 
different sign and is obtained as long as the eigenvalues are in rational independence 
and do not fulfil the condition of "resonance" 

Xm-\i = Q (24) 

for integer vectors m, with \m\ < Umax- 

Analogously, the energy functional (19) is also Taylor expanded and transformed 
according to the change of variables in the equations of motion which results in the 
expansion 

m j 

after introducing canonical coordinates = X2i-i and momenta pi = X2i- Moreover, 
with the latter definition the equations of motion (22)-(23) can easily be integrated to 
a common Hamilton function 



H 



m2i 



™2j 



(26) 
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according to Hamilton's equations if the coefficients j3m satisfy the conditions of 
integrability 

/5m(2i-l) = — /3m(2i); (27) 
Pm{2i-1) (3m{2i'-l) 

= , [Zii) 

m2i m2i' 

l^m{2i) ^ f^m{2i') ^^g) 

m2i-i m2i'-i 

for all i,i' = l,...,d {i ^ i'). In order to guarantee both, the satisfaction of these 
conditions of integrabihty as well as the equivalence of the integrated Hamiltonian H 



with the energy functional (25) an additional transformation is necessary. For this 
purpose, we scale the phase space variables with time-independent functions z/g-(qr,p) 
and i^pi{q,p) according to 

Qi = ^q, {q, p) Qt , Pi = ^p, {q, p) Pi- (30) 
with the constraint of their product 

fii{q,p) = i^gXq,P) ^pAq^P) = 'l + ^f^m YliQjPjr' (31) 

m j 

to be a formal power series of the products qjPj, as well as an appropriate choice of the 
Hm |7] finally guarantees the equivalence of the integrated Hamiltonian (26) with the 
transformed energy functional: 

H (J) = E{J) = H (J) (32) 

Here, action variables J, = qiPi and Jj = iqiPi, respectively, have been introduced 
depending on whether the corresponding eigenvalue of the linearised equations of motion 



is real or purely imaginary. Equation (32) finally serves as classical Hamilton function 
in the sense that it locally reproduces the energy of the system and its Hamiltonian 
equations of motion describe the dynamics in the vicinity of the fixed point equivalently 



to equation (16) 



2.4. Thermal decay rates 

Within the variational approach to monopolar BECs using coupled Gaussian wave 
functions, the set of differential equations (16) exhibits two fixed points [5]. One of 



them is stable corresponding to the metastable ground state of the condensate and the 
other one is of saddle-centre-. . . -centre type corresponding to an unstable excited state. 

These properties, of course, also hold after having applied the near-identity 
transformation described above in order to locally map the variational parameters to 
phase space. The constructed Hamilton function in phase space, thus, takes the form 
depicted in figure [l| featuring a local minimum and a saddle. The latter has precisely 
one unstable direction and can, therefore, be used to divide the phase space into a region 
of "reactants" formed by the metastable BEC and a region of "products" in the form 
of the collapsed condensate. Calculating the decay rate is, thus, possible by applying 
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Figure 1. Schematic drawing of the phase space structure of the constructed Hamilton 
function in equation (32). The metastable ground state of the BEC corresponds to 
a local minimum, and classical decay is possible after thermal excitation. If the only 
decay channel (solid arrow) requires crossing a saddle point in phase space the thermal 
decay rate is given by the Boltzmann average of the flux over this saddle. 



TST (see reference [6]) in association with the constructed Hamiltonian, because the 
only possibihty of the BEC to collapse is by crossing this saddle, and the decay rate is 
given by the flux over it. 

At a fixed energy, the directional flux through the dividing surface between 
"reactants" and "products" is given by (6)|15|[16| 



f{E) = i2ny-'ViE) (33) 

with V{E) being the phase space volume of the actions {J2, ■ ■ ■ , Jd) which is enclosed by 
the contour H{0, J2, ■ ■ ■ , Jd) ^ E and Ji = corresponding to the "unstable direction" 
of the saddle. If the condensate is in contact to a bath of finite temperature the thermal 



decay rate is then given by the Boltzmann average of equation (33). After a short 
calculation, this yields (cf. reference [II] ) 

where /5 = l/k-BT, and Zq is the canonical partition function. Because nearly all states 
will be localized in the vicinity of the ground state, we can well approximate the latter 
by 

Zo = l|dJ(...dJ^e-^^'(^-'^^) (35) 
with H'{J[, . . . , J'^) being the normal form expansion at the local minimum. Furthermore 



taking into account the particle number scaling (10), the thermal decay rate is given by 



1 [e-^W,^2,...,J.)dJ2...dJ^ ^ ^ 

t, (36) 



27r/3/dJ( ...dJ^ e-^W-^i' 
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where we identify as the particle number scaled inverse temperature. 



However, both integrals in equation (36) will, in general, not converge, which is due 



to the fact that the normal form expansion has been truncated at the order nmax- We 
will therefore restrict the area of integration to the condition 
dH 

^. = > (37) 

for all i, in view of the fact that all frequencies occurring on the tori in phase space have 
to be non-negative. 

3. Results and discussion 

In order to calculate the thermal decay rate of BECs with monopolar interaction, we 
first determine the stable and the unstable fixed point of the equations of motion. 



equation (16), for given physical parameters iV 7 and Osc and Taylor expand these in 



the vicinity of the fixed points up to a chosen order nmax- Then, we proceed as described 



in section 2^ to map the variational parameters to classical phase space variables and 
obtain the corresponding local Hamilton functions II{J) and H'{J'), respectively, from 
the transformed mean-field energy functional. With their knowledge, the decay rate is 



calculated from equation (36) under the constraint (37) 



As has been shown by Rau et al |5] the main contribution from the extended 
variational ansatz occurs when the number of wave functions is increased from = 1 to 
Ng = 2. We will therefore restrict ourselves to that case in the following and, moreover, 
to Hamiltonians up to fourth order of the action- variables where we observe convergence. 

For Ng = 2 coupled Gaussians we have eight real variational parameters x, one 
of which is fixed by normalizing the wave function to J d!^r \tlj{r,t)f = 1 and another 
corresponds to a global phase that can be set to zero, so that we are left with six 
independent ones. Determining the corresponding classical Hamiltonian in fourth order 
approximation in J is already non-trivial since, in this case, the expansion of the mean- 
field energy functional up to eighth order in the variational parameters x (scalar valued 
polynomial with 3003 terms) and that of the equations of motions in seventh order of x 
(vector valued polynomial with 10 290 terms) are required. After the mapping to phase 
space, these are simplified to a fourth order polynomial of J with 35 terms, which is a 
reduction of the number of monomials by altogether 99.74%. 

Note that a monopolar BEC, as investigated in this paper, features the phenomenon 



of self-trapping under certain conditions 10 , i.e. an external trap is not necessary to 



keep the condensate stable. However, at least a weak trap is required here to avoid a 
dissolving of the BEC, which otherwise would be a second decay channel. An external 
trap, thus, guarantees that the only decay mechanism of the BEC is its collapse. 

Figure |2] shows the thermal decay rate of a monopolar BEC calculated in third 
order normal form in J using a single Gaussian trial wave function (solid line) and 
Ng = 2 coupled Gaussians (dashed-dotted line). For the calculation, we used a particle 
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Figure 2. Comparison of the thermal decay rate of a monopolar BEC described 
with Ng = I (solid line) and Ng — 2 (dashed-dotted line) Gaussian wave functions 
in dependence of the scattering length Ogc in third order normal form of the action 
variable. The data have been calculated for a particle number scaled inverse 
temperature of N'^(3 — 900 and a trap frequency of N^'y'^ ~ 1.0 x 10^'^. It can be 
seen that the critical scattering length is shifted to higher values when increasing the 
number of Gaussians and that, for a fixed scattering length, the decay rate rises by 
several orders of magnitude. 



number scaled inverse temperature of = 900 and a weak trap with a frequency of 

One consequence of the use of coupled Gaussian wave functions is that the critical 
scattering length Ocrit below which the condensate cannot exist is shifted to larger values 
(cf. reference |5|). For the parameters used here, this is the case from a^nt ~ —1.145 for 
a single Gaussian trial wave function to Ocrit ~ —1.024 for the two coupled ones. Figure 
[2] reveals this behaviour for the whole curve and shows that the thermal decay rate 
calculated with a single Gaussian, as it has also been used in references [2||3] for BECs 
without long-range interaction, underestimates the result of the extended variational 
ansatz by several orders of magnitude for a fixed value of the scattering length age- 
Considering the point of the critical scattering length, the decay rate changes only very 
little compared to a single Gaussian, and the general dependence of the decay rate on 
the scattering length is retained exhibiting a rapid monotonic increase when decreasing 
the scattering length. This increase, however, becomes weaker when one approaches the 
critical value. 

Figure |3^ shows the thermal decay rates for different normal form orders. The 
calculations have been performed for Ng = 2 coupled Gaussians and for the same physical 
parameters used in figure |2} The first-order approximation (dashed line) overestimates 
the decay rate over the whole range of the scattering length, whereas using the normal 
form Hamiltonian in second order in J (dashed-dotted line), we observe the smallest 




-0.5 I ■ ' ■ ' ■ ' ■ ' ■ ' ' ' — 

-1.04 -1.02 -1.00 -0.98 -0.96 -0.94 -0.92 -0.90 

scattering length a 

Figure 3. a) Thermal decay rate of a nionopolar BEC described with Ng = 2 coupled 
Gaussian wave functions in dependence of the scattering length Osc and normal form 
orders (NFO) 1 to 4 of the action variables. Temperature and trap frequency are the 
same as in figure [2] Right inset: The thermal decay rates obtained from the third- and 
fourth-order normal form Hamiltonian cannot be distinguished any more, indicating 
convergence. Left inset: At Cgc ~ —0.999 the eigenvalues are close to "resonance" 
and the normal form expansion as well as the decay rate diverge. As shown in b) 



this is the case because the condition of resonance (24 1 is numerically fulfilled for 
m ~ (0, 0, 1, 0, 6, 0) for this set of physical parameters. 



values throughout. However, the results calculated by the third- and fourth-order 
Hamiltonian (dotted and solid lines) cannot be distinguished within the line width of 
the plot (right inset in figure [3^), indicating convergence. 

At a scattering length of about Osc ~ —0.999, one observes a strong deviation of 
the calculated decay rate in the fourth order approximation from all the other curves 
(left inset in figure [3^), which is in contrast to the behaviour all along the rest of the 
investigated range of the scattering length. As shown in figure [Sjo, the eigenvalues Aj 
of the linearised equations of motion which are used for the normal form expansion 
run into "resonance", i.e. equation (24) is fulfilled within the numerical accuracy for 
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Figure 4. a) Thermal decay rate of a monopolar BEC for a fixed scattering length 
of flsc ~ —0.96, a trap frequency of N'^'j'^ = 1.0 x 10"'^ and normal form orders 1 to 4 
in dependence of the particle number scaled temperature N'^p. The BEC is described 
with Ng ~ 2 coupled Gaussian wave functions, b) The relative deviation Si defined in 



equation (38) is used in order to estimate the convergence of the procedure. 



the integer vector m = (0,0,1,0,6,0) in seventh order of the variational parameters, 
\m\ = 7 (corresponding to the fourth order in J after integration). This leads to the 
divergence of the fourth order normal form Hamiltonian and with it the decay rate at 
Osc ~ -0.999. 

Moreover, the convergence behaviour of the decay rate with increasing normal 
form order strongly depends on the temperature of the system (see figure |4^). While, 
we observe fast converging results for low temperatures and large particle numbers, 
respectively, i.e. large values N'^P, where the decay rates calculated from the third- and 
fourth-order normal forms match within the line width for N'^P > 800 the convergence 
becomes worse when decreasing the scaled inverse temperature N'^13. For N^/S < 200 
the calculations even show a monotonic increase of the decay rate with higher normal 
form order. 

In order to estimate the convergence of our results, we use the relative deviation 

= (rNFO=i ~ rNFO=i-l)/rNFO=j-l (38) 

shown in figure [Ija. The corrections to the decay rate obtained from the second- (i = 2) 
and the third- (i = 3) order normal form are significant throughout, while this is only 
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true for low N'^(5 for the fourth order (z = 4). For large values of N'^(5 the corrections 
quickly shrink and in case of > 1000 these are of the relative order of 10~^ to 10~^, 
clear evidence again of the convergence of the decay rate in fourth-order approximation 
in the low-temperature regime. 



4. Summary and outlook 

We have demonstrated the applicability of a variational approach to classical transition 
state theory by means of calculating thermal decay rates of Bose-Einstein condensates 
with an additional long-range interaction using coupled Gaussian wave functions. 
This procedure has proven as a powerful tool for that purpose: For the extended 
variational ansatz with coupled Gaussians, we observed convergence in eighth order of 
the variational parameters at cold temperatures and high particle numbers, respectively. 
Moreover, the results show that previous estimations using single Gaussian wave 
functions [2]|3], on the one hand, reveal a good qualitative agreement with those obtained 
from the extended ansatz but, on the other, underestimate the thermal decay rate by 
several orders of magnitude for a fixed value of the scattering length. 

To further improve the results and to also achieve convergence for higher 
temperatures the procedure can be extended to the use of more than two coupled 
Gaussian wave functions and higher normal form orders. This is also necessary in 



order to compute decay rates of experimentally accessible dipolar BECs 17 where the 



symmetry breaking dipole-dipole interaction and the occurrence of blood-cell shaped 



condensates 18 require up to six coupled and non-radially symmetrical Gaussians to 
reach the accuracy of the numerical results [5]. Convergence of the decay rate is, finally, 
expected when increasing both the number of Gaussians as well as the normal form 
order. 
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